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Abstract 

This paper generalises King's models for tidally truncated star clusters by including 
approximately the non-spherical symmetry of the tidal field and the resulting non-spherical 
distortion of the cluster. 
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1. Introduction 

The problem of constructing a model of a stellar system such as a globular cluster has 
a long history. In the last twenty years or so, however, much the most widely used models, 
at least in this context, are King's models (King 1966) and various extensions of them. 

Several factors have led to the popularity of these models. First, they are very easy 
to compute. Second, they are successful. Even the single mass models are an excellent 
fit to the light distribution of most observed galactic globular clusters (e.g. Illingworth & 
Illingworth 1976). Third, they are easily extended (with suitable assumptions) to include a 
spectrum of masses or anisotropy in the distribution of velocities, and indeed such models 
are the standard means of interpreting radial velocity data (e.g. Gunn & Griffin 1979, 
Meylan 1988). Fourth, they approximately incorporate essential aspects of theoretical 
stellar dynamics: (i) Jeans' theorem (ii) relaxation effects and (iii) tidal truncation. 

It is the aim of this paper to remove some approximations involved in these models, 
especially with regard to (i) and (iii). In the construction of King's models the (specific) 
energy of a star is taken as v 2 /2 + $i where v is the speed and $i is the potential due 
to the cluster. In fact, each star is moving in the potential field of both the cluster and 
the galaxy. Even so, the energy is not conserved, because v is measured in an accelerating 
frame and centrifugal terms must be added. 

In this paper, then, we shall construct a series of models analogous to King's series, 
but approximately making allowance for the non-spherical symmetry of the tidal field and 
the non-conservation of energy. The method we use is an expansion of the potential in 
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spherical harmonics, though truncated at very low order. We also use the same 'lowered' 
Maxwellian distribution function as in King's models, except that it is now a function of 
the Jacobi integral, which is a conserved quantity when the potential is in uniform rotation. 
Thus, we envisage a cluster whose bary centre rotates on a circular orbit in the plane of 
symmetry of an axisymmetric galaxy. A shortened version of our results is presented in 
Heggie & Ramamani 1993, and we would like to mention that a similar set of models is 
discussed by Weinberg (1993). 

The outline of this paper is as follows. In section 2, we present the necessary prelim- 
inaries. Section 3 describes the numerical method. In section 4 results are discussed for 
two cases (i) when the galaxy is modelled as a point mass and (ii) when the tidal field is 
modelled using observational data from the vicinity of the sun. 

2. Theory of a tidally perturbed model star cluster 

2.1. Case of a point mass perturber 

The cluster under consideration has mass Mi, core radius r c and a limiting radius of 
order r e . It contains stars of equal mass and these stars will have an isotropic velocity 
distribution with respect to a rotating frame of reference which will be described shortly. 
The cluster moves in a circular orbit of radius D around the Galaxy, which is represented 
by a point mass of mass M 2 . The changes in the Galactic potential across the cluster 
are small for D ^> r e , and the field of the Galaxy is adequately represented in a linear 
approximation. 
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We choose a coordinate system rotating with an angular velocity equal to that of the 
orbital motion of the Galaxy-cluster system, and choose its origin to be at the centre of 
mass of this system. The x-axis is chosen to be pointing towards the Galaxy, the y-axis 
perpendicular to the x-axis in the orbital plane in the sense of the orbital motion, and 
the z-axis perpendicular to the orbital plane in the right handed sense. We denote all 
attributes of the cluster and the Galaxy with the suffices 1 and 2, respectively, and the 
central values of the cluster with the suffix c. 

A star in the cluster, at a point r = (x,y,z), has an integral of motion J, the Jacobi 
constant, which is given by 

J = $ 1 + $ 2 _Io;V + y 2 ) + ^ 2 , (1) 

where <E»i, $ 2 are the potentials at r due to the cluster and the Galaxy, respectively, uj is 
the angular velocity of the frame of reference, and v is the speed of the star in the rotating 
frame. Writing the kinetic energy per unit mass as T = \v 2 , we have J = U + T, where 

U = <f> 1 + <f> 2 -±w 2 (x 2 + y 2 ) (2) 

- *i - ^ - GM 2 jD\x\ - \y\ - \z 2 ) - ^ - GM 2 /2D^x\ + y 2 ) (3) 

where xi,yi,zi are the coordinates of the star with respect to the centre of mass of the 
cluster, with the directions of the coordinate axes chosen as before. Expression (3) is 
obtained on the assumption that M 2 ^> Mi, and r\ <^ D. The constant terms here are 
irrelevant in what follows, by redefinition of the zero-point of the potential. 

From Jeans' theorem, the phase-space distribution function (mass- density) for the 
stars of the cluster can be taken as any function of J, and in order to make as close an 
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analogy as possible with King's models we take 

/(r, v) = K(exp(—2j 2 J) - exp(-2j 2 J e )), (4) 

where J e is the value of J at the limiting zero velocity surface, at the edge of the cluster, 
and K, j are constants. Then the density distribution is, as usual, 

p= [ * f(r,v)4jrv 2 dv, (5) 
Jo 

where v e is the maximum speed at the point r = (x\, y\, zi), which is given by 

v 2 e =-2(U-J e ). (6) 



We now exploit the arbitrariness in the zero-point of the potential by supposing that 
J e = 0. Then it is easy to see that the density takes the form 

p = p c R(j 2 U) : (7) 

where p c is the central density and R is expressible in terms of the error function, exactly 
as in the usual theory of King models (King 1966). Where the theory differs is in the 
expression for U, which may be written 

U = $ l +Q 1 (8) 

where Q is the sum of the quadratic terms in expression (3). Now if these terms were 
neglected we would obtain King's models, and $i would be the corresponding potential. 
If the tidal field, represented by Q, were weak, then $i would differ slightly. We therefore 
write 

$i = $ fc + 5$, (9) 
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where $fc is the potential of the corresponding King model, and 5$ is the change in the 
potential due to the distortion of the cluster by the tidal field. As usual, potential and 
density are related by Poisson's equation: 

V 2 $i = 4nGp(U). (10) 

Here is another difference between these models and King's: p depends on U, not on <£>!. 

Now we take the important step of linearising with respect to the deviation from the 
King model, writing 

p(U) = p($ k ) + (5$ + Q)-^- + .... (11) 
Then from eqs. (9) - (11) it follows that 

V 2 d> fc = 4-kGp^u) (12) 

and 

V 2 M> = 4nG^P±{5$ + Q}. (13) 
Eq. (12) is the same as eq.(12) in King (1966). We now write 



5<S>(r,e, ( f>) = J2^n(r)S n (6, ( f>) (15) 



and 

Q = Y J Qn{r)S n {9A), (16) 

where r, 6, 4> are spherical polar coordinates with origin at the cluster centre, and the 
S n (9, 4>) are a complete set of surface spherical harmonics. Eq.(13) can then be written as 

[V 2 - F{r)] S$n(r)S n (e, <f>) = £ F{r)Q n {r)S n {9, </>), 
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which implies 



V 2 M> n (r) - n(n r + 1} ^ n (r) - F(r)M> n (r) = F(r)Q n (r). (17) 



In fact only a few terms of these series are needed. We have 

x\ = r sin cos 0, y\ = r sin 9 sin 

and 



zi = r cos 0, 

whence eqs.(3) and (8) give 



Q = -^r 2 { l -S {e^) + S 2 (0,0)}, (18) 



where So (0,0) = 1 and 



S 2 (6,4>) = -^P 2 °(cos0) + ^P 2 2 (cos0)cos20. (19) 
Here P™ are Legendre's associated polynomials of degree n and order m. 



Since Q has only zero- and second-order spherical harmonics, we get only two relevant 
equations, with n = and n = 2, from eq.(17). They are 

® + ?^. 4 A = . 4lC MM«i! ( 20 ) 

dr 2 r dr d$ fc d§ k D 3 3 

and 

« + ^-{«H lG m» I = - 4lG ^).^. (2!) 
dr 2 r dr r 2 d§ k d§ k D 3 

Since So(9,4>) = 1 the quantity 5<E>o is the spherically symmetric part of the change in 

the potential which arises due to the distortion of the cluster, and 5<E>2 multiplied by S 2 



gives the non-spherical part. It is worth pausing here to note the shape of this distortion. 
Indeed from eq.(19) we have 

5 3 5 
S 2 (0 6) = — + - sin 2 6 cos 26- - cos 2 6. 
12 4 4 

Thus S 2 is largest along the xi-axis and smallest along the zi-axis. The relevant values 
are S 2 (n/2,0) = 7/6, S 2 (n/2, tt/2) = -1/3 and 5 2 (0,0) = -5/6. 

2.2 Jacobi integral for a galactic disk 

Chandrasekhar (1942, p. 217) gives the equations 

m(i - 2u c f] + aiO = -dQ/d£, 

771(77 + 2cu c £) = —dfl/dr), 

m(( + a 3 () = -dn/d(, 

for the motion of a star of mass m moving in the potential Q of a cluster. The coordinates 
(^,?7,C) have the same meaning as (xi,yi,2i) in §2.1, cu c is cu, «i and 03 are constants, 
and subscripts i have been dropped. If O = 0(^,77, £), i.e. we suppose the potential due 
to the cluster is time-independent in the rotating frame, then multiplying these equations 
by £, ?7, C, respectively, and adding, gives the Jacobi integral 

J = Id 2 + V 2 + C 2 ) + + \a,e + n/m. (22) 

Denoting by B the galactic potential, Chandrasekhar (p. 217 again) gives the definitions 

a 1 = ^E-l^ (23) 
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«3 = 0, (24) 

where z is height above the galactic plane and zu is the distance from the symmetry axis 
of the galaxy. Using Poisson's equation in cylindrical polars (eq. 4.528 in Chandrasekhar) 
we find that 

d 2 B 1 dB 

a 3 = AixGpc - j— — , 

ow z zu ozu 

where p G is the galactic density. From the first of eqs. (5.501) in Chandrasekhar we deduce 
that 

26 dO 



as = AjiGpo 



zu dm ' 

where 6 is the local speed of circular galactic motion. Hence, by eq.( 1.429) in Chan- 
drasekhar, 

0:3 = AirGpa — 2lu(A + B), (25) 
where A and B are Oort's constants. Similarly 

ai =4A(B - A), (26) 

by Chandrasekhar's eq. (5.615). 

The Jacobi integral has exactly the same form as in the case of a point mass Mq at 
distance D. In this case B = -GM g /(vj 2 + z 2 ) 1/2 and p G = 0, whence a-y = -3GM G /D 3 
(by eq.(23), where the potential derivatives are evaluated at z = 0, w = D) and CI3 = 
GMq/D 3 . Thus the Jacobi integral is 

3 = 2 (C + 71 + C ) ~ 2^ + 2^ C + ° /m - 
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This is essentially the same result as eq.(3) of section 2.1. 

Returning to the general case (eq.(22)) and defining Q by analogy with eq.(8), we 

have 

where £, ( have the same meaning as x\, z\ in §2.1. Hence 

Q = lc n r 2 (^S + S 2 ) 

where now 

5o = |(l + ^)P ° (27) 
2 cti 

and 

5 2 = (^-i)P 2 ° + ip 2 2 cos2^. 
ai 2 4 

3. Numerical method 

In order to construct models it is necessary to solve eqs.(12), (20) and (21) simulta- 
neously throughout the cluster by numerical methods. Following King (1966), we choose 
scaled variables W = —2j 2 &k and R = r/r c , where r c is the core radius defined be- 
low. Analogously the perturbations in the potential are scaled to give ys = — 2j 2 5$ and 
y 5 = — 2j 2 5$2- For the sake of definiteness we consider the case of §2.1, and define 

M 2 

Pi - 4 n3 ' 

which is a measure of the density of the perturbing Galaxy. Following King again, we 
take the core radius as satisfying the relation 87rGj 2 p c r 2 = 9. With R as the independent 
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variable, the equations then give 



and 



d 2 W 2 dW _ p(W) 

dR 2 R dR ~ 9 p c ' 1 8J 

d 2 y 3 2 dy 3 n d(p/ p c ) 2 pi d(p/ p c ) 

~dR 2 R~~dR 9 ^M^ V3 = ~ 9R J C -^T' (29) 



d 2 y 5 2 ^ _ 6 + 9 <M^^ = _^ r2Pi d_{p]j>l 



dR 2 R dR R 2 * dW yo p c dW ' 

where, following King, 



p exp(VF) J Q exp(-r])r] 3/2 dr] 
Pc exp(W c ) J^ c exp(— rj)rj 3 / 2 dr] 



(31) 



and 

d jp/Pc) = P_ + 

dW p c eMW c ) f™ c exp(-r])r] s / 2 dr] 



\f _ U \P Pc P, ^_ /o 9 N 

(p/pc) = —^7- = — + — rw c — ; — ; ~ • ( 32 ) 



Here, as usual, a subscript c denotes a core or central value. 

Since p depends only on W, eqs.(28) to (30) may be conveniently transformed with 
the choice of W as the independent variable. We also use ln(l + R 2 ) rather than R as 
the dependent variable in transforming eq.(28), because this varies nearly linearly with W 
both near R = and in an isothermal halo. Thus the equations to be integrated take the 
form 



dyi 
dW 

dy2 
dW 

dm 
dW 

dy *- Hl + R2) \l\-R^ d ^-y/-^ + ^yA (36) 



= V2 (33) 
= (H) 2 (6 + 2R 2 + 9^(1 + R 2 fy 2 ^ (34) 
= 2/4 (35) 



dW 4 R 2 * A \ p c dW y ° dW p c 
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w = m and (37) 

dye _ 3 (1 + R 2 ) 2 2 f 2y 5 2 pi d(p/p c ) d(p/p c ) p \ 



where y\ = ln(l + R 2 ) and the additional variables j/2,J/4 and j/6 are essentially defined 
by eqs.(33), (35) and (37). Thus the original set of three second-order equations has been 
transformed to a set of six equations of first order. 

The numerical treatment of these equations is straightforward. The above system can 
be written as 



dy . , 



and can be discretised as 



Yt+1 Ax yt =g(*,y), (39) 



where x = (xi +Xi+i)/2, etc., and a subscript % indicates the value at the i-th mesh point. 

The boundary conditions at the centre of the system can be obtained using the fol- 
lowing developments. Let 



W - W c = aR 2 + bR 4 + . . 



y 3 = a R 2 + b R 4 + . . . (40) 

y 5 = a 2 R 2 + b 2 R 4 + ... 
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Substitution of eqs.(40) into eqs.(28 - 30) is found to give the results 

a = -3/2, 

a = 0, 

b ° = ~™—(p/Pc)\ and 
zu p c 

6 2 = -^(p/ Pc ) / (a 2 + 3^), 

J-^ Pc 

where (p/p c )' is evaluated at the centre, and a 2 is yet to be determined. (It must be 
solved for along with the solution of the differential equations.) From these results it is 
easily found that the boundary conditions at the centre are y 2 = —2/3, y& = — 2a 2 /3, and 
7/1 = 7/3 = 2/4 = 2/5 = 0. At the outside of the cluster 5$ 2 should join smoothly to the 
exterior (vacuum) solution proportional to -R -3 , whence 

3 1 + R 2 

Eqs.(39), along with the above boundary conditions, were solved numerically using the 
Newton-Raphson method. 

As a check of the numerical method, at least as far as the equations for yi and y 2 
are concerned, King's models were regenerated using our code, and good agreement was 
obtained with the results of King (1966). This yields a one-parameter family of solutions, 
depending on W c . But the system of eqs.(33)-(38) also possesses another dimensionless 
parameter, viz. p\j p c . In physical terms the number of parameters can be understood by 
considering the number of radii which characterise each solution. In fact there are three: 
the core radius, r c ; the radius of the system, r e ; and, finally, the radius of the last closed 
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zero-velocity surface, on which the Lagrangian points of the cluster-Galaxy potential lie. 
These three radii give two independent dimensionless ratios, i.e. two parameters. 

A solution of the equations having been obtained for given values of the two parameters 
W c and pi/pc, the functions <E>fc, 5<E>o and 5<E>2 can be determined from the runs of 1/3 and 
2/5 against W, and then $1 can be found from eqs.(9) and (15). The limiting surface is 
obtained by setting U = 0, where, as before, U = $fc + 5& So + + Q- 

Of greatest interest are those systems in which the edge of the cluster coincides with 
the last closed zero-velocity surface. For a fixed value of W c this can be found by iterating 
the value of pi/p c until a point on the xi-axis is found at which both U and dU/dx\ 
vanish. (This point is the Lagrangian point.) 

We now discuss the special features of the problem discussed in §2.2. In the disk model 
of the Galaxy, the values of «i and as were calculated at sun's distance using A = 0.0144 
km/sec/pc and B = —0.012 km/sec/pc with the density at the solar distance in the Galaxy 

p G = 0.11M o /pc 3 , 

(cf. Kuijken & Gilmore 1989). Using eqs.(25) and (26) we obtain 

^ = -3.99, 

and so So = —4.49, by eq.(27). Again S2 is largest along the xi-axis and smallest along 
the zi-axis: the relevant values are ^2(^/2, 0) = 3.00, £2(71-/2, ir/2) = +1.50 and 62(0, 0) = 
—4.49. Comparison with the values of S2 given for the point-mass model indicate that the 
effects will be more pronounced in the case of the disk model. 
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4. Results and Discussion 

4.1 Description of the models 

Space density profiles of the resulting models are shown in Figs. 1-3 and Figs. 6-8 for 
the point-mass and the disk model of the galaxy, respectively. The space density profiles 
along the ^i-axis are compared with the corresponding King models in Figs. 1 and 6. Figs. 
4, 5, 9 and 10 show the surface densities of our model clusters along the x\- and zi-axes 
as seen from the y\ direction. It can be seen from the results that: 

1. The clusters are triaxial, with the longest axis in the direction to the galactic 
centre and the shortest in the direction at right angles to the orbital plane. Denoting by 
a, 6, c the semi diameters in the x±, y\ and z\ directions, respectively, it is found that the 
axial ratios a/b and b/c decrease only slightly as the scaled central potential W c increases 
from 2.5 to 10. Typical values are a/b = 1.5, b/c = 1.04 for the point mass perturber and 
a/b = 1.50, b/c = 1.34 for the tidal field representing a disk galaxy. The strong flattening 
in the z\ direction in the latter case is already known from iV-body simulations (Terlevich 
1987). The axial ratio will clearly depend on the values of 0:3/0:1. 

2. The ratio of the maximum extent of the cluster to that of the corresponding King 
model decreases slightly as W c increases and is almost the same in the case of a point mass 
perturber as in that of a disk galaxy: a representative value is 1.5, i.e. the approximate 
self-consistent model is 50 percent larger. Interestingly, a similar conclusion was reached 
by Spitzer(1987) for a model in which both the cluster and the galaxy are treated as point 
masses. This agreement suggests that departures of the cluster potential from spherical 

16 



symmetry are modest and helps to justify the drastic truncation of the series in equation 
(15) at n = 2. 

3. The central regions of the clusters are almost unaffected; they are nearly spherical 
and differ little from the central regions of the corresponding King models. 

4.2 Discussion 

Our result (2 above) that our models are more extended than the corresponding 
King models may have significant consequences for the application of King models to 
observations of globular clusters. Since limiting radii have sometimes been inferred by 
fitting King models, it is clear that the resulting estimates are likely to be systematically 
too small. Furthermore, since the inferred strength of the tidal field or the inferred mass 
of the cluster (Freeman 1980, Kontizas & Kontizas 1983) depends on the cube of the 
limiting radius, still larger systematic errors are liable to be present in the estimates of 
these quantities. Nevertheless, it should not be forgotten that the assumptions underlying 
both King models and the models discussed in this paper are not strictly applicable to 
globular clusters at all, since their galactic orbits are unlikely to be circular. 

From the theoretical point of view, our models may be useful for the construction of 
initial conditions for iV-body models of tidally perturbed clusters. At present the initial 
conditions used in such simulations are spherically symmetric, and over the first few cross- 
ing times the system presumably settles into a quasi-equilibrium in the non-spherically 
symmetric tide. 
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In order to check these ideas a number of 1000-body simulations have been conducted 
using NB0DY5 (Aarseth 1985) with three different sets of initial conditions and parame- 
ters: (i) one of the models described in this paper, with scaled central potential W c = 2.5, 
and the tidal parameters given at the end of §3; (ii) a King model with the same central 
potential and tidal parameters; and (iii) the same King model without an external tide. 
All models have stars of equal mass, and units are standard (Heggie & Mathieu 1986), 
with the crossing time being 2y/2 units. Figs.(ll)-(13) show, for these three models, the 
mean square values of x\, y\ and zi, measured relative to the density centre (Casertano 
& Hut 1985). 

These results indicate that the model constructed according to the theory of the 
present paper remains in approximate dynamic equilibrium for about the first 4 time 
units. (The dip in one of the curves in Fig. 11 is no greater than that in one of the curves 
in Fig. 13, indicating that it can be ascribed to finite- N effects). After about t ~ 4 there 
begins a gentle expansion (on a time scale of about 70 time units). The tidally perturbed 
King model, on the other hand, flattens in the zi-direction over the same period of about 
4 units and spreads in the other two directions. Indeed, at about t = 4 it quite closely 
resembles the model shown in Fig.ll. Eventually this model shows the same tendency to 
expand, and it does so on a very similar time scale to the model shown in Fig.ll. We 
interpret the expansion as being due to two-body relaxation, and this is confirmed by the 
isolated model, which shows no tendency to triaxiality, but does expand in the last part 
of the run. 

In conclusion we remark that it would be desirable to improve the models by a higher- 
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order expansion, or use of a grid-based potential solver. Such a method has been used by 
Rix and White (1989) to discuss models for a binary galaxy. The models obtained by us 
can be regarded as a limiting and simplified case of the same problem. 
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Figure Captions 



Fig.l Space density on the xi-axis for nine of the models described in this paper and the 
corresponding King models, in the case in which the tidal field is due to a large distant 
point mass. The scaled central potentials W c (denoted W on the figure) of the models are 
2.5, 3, 4, 5, 6, 7, 8, 9 and 10, and the King models are distinguished by dotted curves. 
The scaled distance from the cluster centre is R, and the density is scaled by the central 
value. 

Fig. 2 As Fig.l, but for the yi-axis. The King models are not shown. 
Fig. 3 As Fig.l, but for the zi-axis. 

Fig. 4 The projected (surface) density of the nine models displayed in Figs. 1-3. The surface 
density is scaled to the central value. The cluster is viewed along the yi-axis and the profile 
along the x i -axis is plotted. 

Fig. 5 As Fig. 4, except that the profile along the zi-axis is plotted. 

Fig. 6 As Fig.l, except that the tidal field is a model of that in the solar neighbourhood 
("disk" case). 

Fig. 7 As Fig. 2, but for the disk case. 
Fig. 8 As Fig. 3, but for the disk case. 
Fig. 9 As Fig. 4, but for the disk case. 
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Fig. 10 As Fig. 5, but for the disk case. 

Fig. 11 Evolution of a 1000-body model with initial conditions obtained from one of the 
models described in this paper. The model had scaled central potential W c = 2.5 and the 
tidal field is a model of the tide in the solar neighbourhood. The quantities plotted are 
the mean square values of the coordinates xi, yi and z\, but using as origin the "density 
centre". The lowest curve corresponds to zi, the uppermost to x\. 

Fig. 12 As Fig.ll, except the initial conditions are constructed using a standard (spherically 
symmetric) King model. Again the lowest curve corresponds to z±, the uppermost to x\. 

Fig. 13 As Fig. 12, but without any tide. The dashed curve gives one third of the mean 
square distance from the density centre, (x\ + y\ + z\ )/3. 
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